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I.  Introduction: 


An  important  problem  in  image  processing  applications  1b  the  segmentation 
of  an  image  field  into  disjoint  regions  which  may  possess  the  same  average  gray 
level  but  differ  in  the  spatial  distribution  of  gray  levels.  These  two  character¬ 
istics  are  generally  referred  to  as  tone  and  texture  respectively,  although  more 
precise  definition  of  these  terms  has  remained  elusive.  It  has  been  emphasized 
by  Haralick  [l],  [2],  among  others,  that  a  subtle  relationship  exists  between 
tone  and  texture  which  depends  very  much  upon  the  resolution  with  which  an  image 
is  viewed.  At  both  low  and  high  resolution  the  dominant  feature  is  that  of  tone, 
while  at  intermediate  resolutions  texture  is  often  the  dominant  feature.  The 
most  widely  accepted  definition  of  texture  at  present  [ 3 ]— [ 5 ]  consists  of  a  basic 
local  order  or  quasi-homogeneous  pattern  which  is  repeated  in  a  "nearly  periodic" 
manner  over  some  image  region  which  is  large  relative  to  the  size  of  the  local 
pattern.  We  accept  this  as  a  working  definition  in  what  follows  although  with 
some  qualification  on  the  "nearly  periodic"  repetitiveness  attribute. 

A  number  of  techniques  for  texture  discrimination  have  been  proposed  and 
have  achieved  considerable  success,  although  generally  under  well-defined  and 
rather  limited  operating  conditions.  These  techniques  can  be  classified  as 
either  structural  or  statistical  in  their  approach.  Our  interest  here  will  be 
in  a  purely  statistical  approach.  Structural  approaches  are  described  by 
Zucker  [6],  and  Lu  and  Fix  [7],  among  others.  Recent  work  by  Haralick  [8]  pro¬ 
vides  a  comprehensive  survey  of  most  existing  statistical  techniques  which  he 
classifies  into  eight  broad  categories  possessing  some  degree  of  overlap. 

For  example,  consider  Haralick' s  first  three  categories;  those  based  on  auto¬ 
correlation  functions,  optical  transforms,  and  digital  transforms.  In  reality 


1 


these  techniques  are  all  based  upon  Becond-moment  properties  and  can  be 
collectively  classified  as  such.  The  major  classifications  of  statistical  approaches 
to  texture  classification/discrimination  are  then,  following  Haralick,  those 
based  upon:  second-moment  properties,  edge  density,  spatial  similarity,  spatial 
gray-level  co-occurrence  probabilities,  gray-level  run  lengths,  and  finally 
use  of  tvo-dim8nional  (2-D)  autoregressive  modeling  assumptions.  As  Haralick 
correctly  points  out,  the  spatial  similarity  approach  is  restricted  to  binary 
images  while  the  2-D  autoregressive  linear  estimation  approach  is  severely 
limited  in  the  classes  of  texture  for  which  it  is  useful.  In  particular,  the 
2-D  autoregressive  process  does  not,  except  under  pathological  assumptions, 
exhibit  the  local  pattern  replication  attribute  considered  an  essential  in¬ 
gredient  of  texture.  As  a  result,  the  competing  statistical  approaches  for 
general  application  to  texture  discrimination  are  reduced  to  the  four  remain¬ 
ing  categories  as  enumerated  above. 

Many  of  the  specific  techniques  in  these  remaining  four  categories  are 
based  upon  heuristic  or  ad  hoc  arguments  while  a  comparative  few  have  been 
based  upon  the  rigorous  application  of  statistical  decision  theory  concepts 
under  specific  stochastic  modeling  assumptions.  A  comparative  study  of 
several  of  the  more  frequently  used  statistical  approaches  is  provided  by 
Weszka,  et  al.  I 93 •  In  particular,  the  efficacy  of  various  features,  drawn 
from  these  four  categories,  was  investigated  In  the  context  of  terrain  classif¬ 
ication.  While  results  of  this  nature  are  useful  they  provide  little  guidance 
on  how  the  relative  performance  is  affected  under  various  modeling  assumptions. 

It  is  difficult,  for  example,  to  extrapolate  the  results  of  these  and  similar 
studies  to  applications  beyond  the  specific  data  base  for  which  performance 
has  been  provided.  Clearly,  use  of  more  general  2-D  stochastic  texture  models 


whose  parameters  can  be  easily  related  to  texture  properties  would  remedy 
this  situation. 

Much  of  the  work  on  texture  discrimination  has  been  guided  by  human  visual 
discrimination  studies.  Several  researchers  [10],  [ll]  in  this  area  havfe  conclud¬ 
ed,  although  with  some  qualification  [12],  that  humans  can  effortlessly  differ¬ 
entiate  texture  regions  which  differ  in  second-order  statistics  but  cannot 
discriminate  between  regions  which  differ  only  in  third  and  higher-order  statistics 
This  is  supported  by  the  degree  of  success  achieved  with  texture  discrimination 
algorithms  based  upon  second-order  properties  alone,  such  as  second-moment 
techniques  utilizing  autocorrelation  functions  or  power  spectral,  densities. 
Nevertheless,  for  virtually  all  the  existing  techniques  it  is  possible  to 
contrive  counterexamples  which,  although  effortlessly  discriminated  by  human 
observation,  cannot  be  discriminated  by  algorithmic  means. 

*1* 

In  Fig.  1,  for  example,  we  illustrate  realizations  of  two  random  fields'  which 
are  visually  distinct  yet  possess  identical  autocorrelation  functions  and/or 
power  spectral  densities.  These  fields  cannot  be  distinguished  by  algorithmic 
approaches  based  upon  second-moment  properties  alone.  Similarly,  in  Fig.  2, 
we  illustrate  several  realizations  of  random  fields  which  possess  the  same 
number  of  edges  per  unit  distance  and  are  such  that  the  average  gray-level  run 
length  along  any  randomly  chosen  line  segment  is  identical.  Again  these  texture 
regions  are  easily  discriminated  visually  although  algorithmic  techniques  based 
on  either  edge  density  or  gray-level  run  lengths  alone  cannot  discriminate  the 
various  texture  regions.  Finally,  it  is  possible  to  contrive  random  fields 
whose  joint  probability  density  function  at  two  points  separated  by  a  specified 
distance  d  are  identical,  although  this  need  not  be  true  for  all  values  of  d. 

+  The  parameters  defining  the  2-D  random  fields  in  Fig.'s  1  and  2  will 
be  described  in  a  later  section. 


The  implication  here  is  that,  unless  the  separation  distance  d  is  Judiciously 
chosen,  texture  discrimination  algorithms  based  on  spatial  gray- level  co-occurr¬ 
ence  probabilities  are  incapable  of  distinguishing  visually  distinct  texture 
samples.  Again  this  points  out  the  need  for  texture  discrimination  approaches 
based  upon  rigorous  application  of  statistical  decision  theory  concepts  under 
specific  and  flexibly  parameterized  stochastic  modeling  assumptions. 

In  the  present  paper  ve  describe  a  class  of  2-D  random  fields,  of  which 
the  samples  in  Fig.'s  1  and  2  are  selected  realizations,  which  we  feel  provides 
a  realistic  and  conveniently  parameterized  model  of  texture  in  images.  This 
class  of  random  fields  bear  some  relationship  to  the  random  mosaic  models  for 
texture  described  by  Schachter,  et  al.  [13].  Based  upon  this  stochastic  model 
we  propose  a  new  approach  to  texture  discrimination  which  is  an  approximation 
to  the  statistically  optimum  maximum  likelihood  classifier.  This  approach,  for 
reasons  to  be  described,  makes  use  of  the  spatial  gray-level  co-occurrence  matrix 
introduced  by  Haralick  [l],  [2].  However,  unlike  the  Haralick  approach,  we  do 
not  make  use  of  heuristically  defined  features  for  extracting  texture  information 
from  the  spatial  gray-level  co-occurrence  matrix.  Rather  our  approach  is  based 
upon  a  maximum  likelihood  hypothesis  test  of  the  gray-level  co-occurrence  matrix. 
This  leads  to  a  rather  simple  implementation  as  a  2-D  digital  filtering  operation 
on  the  original  image.  Results  indicate  a  substantial  performance  improvement 
over  competing  approaches. 

After  some  preliminary  comments  on  2-D  random  fields  in  Section  II,  the  con¬ 
struction  and  properties  of  the  stochastic  texture  model  are  described  in  Sections 
III  and  IV  respectively.  The  structure  of  the  maximum  likelihood  texture  dis¬ 
criminator  is  provided  in  Section  V  while  the  approximate  digital  implementation 
is  described  in  Section  VI.  Typical  results  are  illustrated  in  Section  VII  while 
Section  VIII  provides  a  summary  and  conclusions 
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II.  Preliminary  Discussion: 
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We  consider  an  image  as  a  family  of  random  variables  {fx(u>),  xeR  }  ,  or 
a  random  field,  defined  on  some  fixed  but  unspecified  probability  space  (ft,A,P). 
For  convenience  we  suppress  the  functional  dependence  upon  the  underlying  prob¬ 
ability  space  and  consistently  write  f(x)  for  f^Cw).  The  covariance  function 
of  the  random  field  then  becomes^ 

=  E{f(x)f(j£)}  ;  x.^eR2  ,  (1) 

2 

where  E{#}  represents  the  expectation  operator.  If  a  random  field  {f(x),  xeR  } 
possesses  a  covariance  function  invariant  under  all  Euclidean  motions  it  will 
be  called  homogeneous  and  isotropic  (cf.  [lU]  for  definitions).  In  this  case 
the  covariance  function  of  the  field  evaluated  at  two  points  can  depend  only 
upon  the  Euclidean  distance  between  these  two  points  so  that 

E{f(x+u)f(x)}=  Rff( | |u| | )  ,  (2) 

T  2 

where  u  =(u^,u2)  is  an  element  of  R  and  | |u| [  represents  the  ordinary  Euclidean 
norm  defined  in  terms  of  an  inner  product  <*,*>  according  to 

I |u| | 2  =  <u,  u>  *  u*  +  uf  (3) 

By  construction,  the  2-D  random  fields  to  be  described  here  are  of  this 
category.  Furthermore,  they  have  been  explicitly  constructed  so  that  the 
Joint  probability  density  function  (p.d.f. )  of  the  field  evaluated  at  two 
points  likewise  depends  only  upon  the  Euclidean  distance  between  these  points. 
More  specifically,  define  the  random  variables  f^fCx)  and  f^=f (x+u) .  The  joint 
p.d.f.  associated  with  these  two  random  variables,  parameterized  by  the  spatial 
coordinates,  then  satisfies 

ptf^fg-,  x >x+u^  =  p{f1,f2;  ||u||)  ,  (*+) 

■x - 

We  assume  the  field  is  of  second  order  (i.e.,  variances  exist)  and  possesses 
zero  mean. 
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which  is  the  2-D  concept  of  stationarity  [15]  or  invariance  which  will  be  most 
useful  for  our  purposes. 

The  corresponding  power  spectral  density  function  is  given  by 


sff<2>  ■  |, 


RffU  |u|  I)exp{-j  <u,  u  >}du 


(5) 


T 

where  u)  *  represents  a  2-D  spatial  frequency  vector  and  du  is  the 

2 

differential  volume  element  in  R  .  This  expression  can  be  evaluated  up  to 
functional  form  with  the  aid  of  a  theorem  of  Bochner  [16]  with  the  result 


Sff(u)  =  S(fl)  =  2rr  j  1  Rff(A)  JQ(Afi)dA  ,  (6) 

0 

where  ft  =  |  |uj|  ^(uj^+u)*)  ^  represents  radial  frequency.  Here  JQ(  • )  denotes 
the  ordinary  Bessel  function  of  the  first  kind  of  order  zero.  The  quantities 
Sff(»)  and  Rff(*)  are  then  related  through  a  Hankel  transform  [17],  [18] . 

An  important  aspect  of  the  approach  to  texture  discrimination  described 
here  is  the  use  of  a  stochastic  texture  model  whose  second-order  statistics 
are  invariant  under  both  translation  and  rotation.  Various  stochastic  texture 
models  proposed  previously  do  not  possess  this  property.  For  example,  much 
use  has  been  made  of  2-D  autoregressive  models  [19]-[21]  for  texture,  often 
under  a  separability  assumption  in  the  two  orthogonal  spatial  directions. 

These  processes  cannot  possess  second-order  statistics  invariant  under  all 
rigid  body  motions.  On  the  other  hand,  we  feel  strongly  that  a  reasonable 
texture  model  should  possess  this  property.  Texture  should  retain  its  ident¬ 
ity  regardless  of  the  orientation  or  perspective  in  which  it  is  presented. 

Use  of  naturally  occurring  textures  exhibiting  obvious  directional  properties, 
such  as  those  in  the  book  by  Brodatz  [22],  serves  only  to  obscure  this  issue. 
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Had  many  of  these  samples  been  presented  in  a  different  orientation  would  a 


different  texture  category  have  resulted?  We  feel  rather  that  these  examples 
should  be  more  properly  considered  sample  fields  or  realizations  of  a 
random  field  as  defined  here.  While  particular  realizations  may  well  exhibit 
directional  characteristics,  the  ensemble  properties,  at  least  up  to  second- 
order  statistics,  should  be  invariant  under  rigid  body  motions. 

Another  criticism  of  existing  stochastic  texture  models,  such  as  auto¬ 
regressive  processes,  is  the  inability  to  account  for  the  predominant  and 
pronounced  edge  structure  present  in  real-world  imagery.  Often  this  edge 
structure  provides  an  important  aspect  of  the.  rather  imprecise  concept  of 
perceived  texture.  Finally,  and  most  importantly,  existing  stochastic 
texture  models  do  not  provide  the  basic  repetition  of  a  local  order  or  pattern 
generally  considered  [3],  [U]  an  important  ingredient  of  texture.  The  stochas¬ 
tic  texture  models  described  in  the  next  section  remove  many  of  these  object¬ 
ions.  Furthermore,  the  mathematical  tractability  associated  with  these  models 
allows  straightforward  development  of  statistically  optimum  texture  discrim¬ 
ination  algorithms. 

Ill  Construction  of  Random  Field  Models  of  Texture: 

The  class  of  random  fields  to  be  used  as  stochastic  texture  models  can 
be  described  as  marked  point  processes  [23]  evolving  according  to  a  spatial 
parameter.  According  to  this  model  the  plane  is  randomly  partitioned  into 
a  number  of  disjoint  geometric  regions  by  an  appropriately  defined  field  of 
random  lines  which  form  the  boundaries  of  these  regions.  The  density  of 
these  random  lines,  or  edges,  is  defined  in  terms  of  a  rate  parameter  X. 

Gray  levels  are  then  assigned  within  elementary  regions  to  possess  correl¬ 
ation  coefficient  p  with  gray  levels  in  contiguous  regions.  We  describe 
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several  schemes  for  partitioning  the  plane  into  elementary  geometrical 
regions.  Given  a  particular  partitioning  scheme,  the  random  fields  are 
completely  defined  in  terms  of  the  two  parameters  X  and  p.  The  parameter 
X  represents  the  "edge  business"  associated  with  an  image  while  p  is  indic¬ 
ative,  at  least  on  an  ensemble  basis,  of  the  "edge  contrast".  For  p  large 
(in  magnitude)  and  negative  there  is  an  abrupt  almost  black-to-vhite  or 
white-to-black  transition  across  an  edge  boundary.  If  p>0,  on  the  other  hand, 
the  transition  across  an  edge  boundary  is  much  more  gradual.  It  is  relative¬ 
ly  easy  to  define  these  parameters  for  wide  classes  of  imagery  data. 

In  the  present  section  we  describe  the  construction  of  this  class  of  2-D 
random  fields.  Relevent  second-order  properties  are  discussed  in  the  next 
section.  We  begin  with  the  case  where  the  plane  is  partitioned  into  random 
rectangular  regions. 

RfccAuiflo&tft-  PeumtiLtion&t  A  fundamental  role  in  the  construction  of  this  class 
of  processes  will  be  played  by  the  integer-valued  random  field  (N(x) ,x>0> 
which  provides  a  2-D  generalization  of  a  counting  process  [24].  In  particular, 
suppose  the  vector  x  is  obtained  from  x  according  to  x=A  x  where  A  is  the 
unitary  matrix 


A= 


cos  8 
-sin8 


sin8 

COS0 


(T) 


defined  for  some  0e[-ir,Tr].  This  transformation  results  in  a  rotation  of  the 
Cartesian  coordinate  axes  (x^,Xg)  by  0  radians  as  illustrated  in  Fig.  3. 
Consider  now  the  integer- valued  random  field  defined  by 

N(x)  *  K^x^  +  Ng(Xg)  ;  x  >  0  ,  (8) 

ip 

By  the  notation  x  >  0  we  mean  that  x  »(x1,x2)  is  such  that  x^>0,  1*1,2 
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where  9e[-ix,n]  is  chosen  according  to  some  p.d.f.  p(6)  and  (N^(£),  £>0}, 
i*l,2,  are  mutually  independent  1-D  counting  processes.  That  is,  ( £ ) 
represents  the  number  of  events  which  have  occurred  in  the  interval  [0,£]. 

We  will  be  particularly  concerned  with  the  case  where  {N^(£),  £>0} ,  i=l,2, 

are  renewal  point  processes  defined  in  terms  of  their  interarrival  distribution. 

The  random  field  {N(x),  x  >  0)  in  (8)  then  assumes  constant  integer  values 
on  non-overlapping  rectangles  whose  sides  are  parallel  to  the  transformed  axes 
)  and  whose  locations  are  determined  by  the  event  times  of  the  corres¬ 
ponding  point  processes  { ( £ ) ,  £>0},  i=l,2  .  Consider  now  the  random  field 
(f(x),  x  >  0)  which  undergoes  transitions  at  the  boundaries  of  these  elementary 
rectangles.  The  gray  level  assumed  throughout  any  elementry  rectangle  is  zero- 

f  O 

mean  Gaussian  with  variance  a4  and  correlation  coefficient  p  with  the  gray 

levels  in  contiguous  rectangles.  More  specifically,  let  X.  represent  the 

i » J 

amplitude  or  gray  level  assumed  by  the  random  field  after  i  transitions  in 
the  X  direction  and  j  transitions  in  the  X_  direction.  The  sequence  {X.  ,} 
is  assumed  generated  recursively  according  to 


X.  t  ,+p  X.  ,  ,-p2X.  +W 

l-l,j  1— J.,j— X  1,J 


;  i.Jli  » 


(9) 


where  |p|<  l,and{W  }  is  a  2-D  sequence  of  independent  and  identically  dis- 

tributed  (i.i.d. )  zero-mean  Gaussian  variates  with  common  variance  o^=a2(l-p2 )2. 

The  initial  values  X^  q,  Xq  k,£  j>  0  are  jointly  distributed  zero-mean 

Gaussian  variates  with  common  variance  a2  and  covariance  properties  chosen 

to  result  in  stationary  conditions.  An  alternative  interpretation  of  the 

sequence  {X  }  is  as  the  output  of  a  separable  2-D  recursive  filter  excited 
i  J 

by  a  white  noise  field.  It  is  easily  seen  that 


t 


For  definiteness  we  assume  Gaussian  statistics.  This  assumption  is  not 
critical  to  the  development  which  follows  and  is  easily  removed. 


Typical  computer-generated  realizations  of  the  resulting  random  field  are 
illustrated  in  Fig.  4  for  selected  values  of  p  when  p(0)  is  uniform  over  [  —it ,tt ] 
and  (Ni(4),  £>0}  ,  i*l,2,  are  Poisson  with  intensities  X^=X,j=X.  The  displayed 
images  here  and  throughout  this  paper  are  square  arrays  consisting  of  256 
elements  or  samples  on  a  side.  In  Fig.  4,  X  is  measured  in  normalized  units 
of  events  per  sample  distance  so  that  there  are  on  average  256X  transitions 
along  each  of  the  orthogonal  axes.  Similarly  in  Fig.  5  ve  illustrate  real¬ 
izations  of  the  resulting  random  field  when  the  point  processes  (N^(£) ,  £>0} , 
i»l , c .  undergo  jumps  of  unit  height  at  equally  spaced  intervals  £=1/X.  The 
starting  positions  ,  i=l,2,  will  be  assumed  uniformly  distributed  over  the 
interval  [0,£j. 

The  preceding  two  examples  are  special  cases  of  the  situation  vhere  the 
point  processes  {N^(£),  £>o),  1*1,2,  are  stationary  renewal  processes  [25],  [26] 

with  Gamma  distributed  interarrival  times.  This  class  of  random  fields  re¬ 
presents  a  2-D  generalization  of  the  class  of  1-D  processes  described  in  [27]. 

In  particular,  ve  assume  the  common  interarrival  distribution  of  the  two 

mutually  independent  point  processes  (N^(£),  £^0),  i*l,2,  possesses  p.d.f. 

v-1 

f(x)  ■  — - —  exp{-x/B)  .  (ll) 

T(v)8V 

vhere  v*l,2,...,  and  0»l/Xv  for  fixed  X>0.  For  example,  if  v*l  ve  have  the 
exponential  distribution 

f(x)  *  Xe  *X  ;  x>0  >  (12) 

associated  vith  the  Poisson  process,  while  in  the  limit  v*»  ve  have 

f(x)  *  6(x-l/X)  ;  x>0  ,  (13) 

corresponding  to  the  case  of  periodic  partitions  as  illustrated  in  Fig.  5* 


In  Fig.  6  we  illustrate  selected  realizations  of  the  resulting  random 
field  for  several  values  of  v  all  with  A=0. 05  and  p=0.0.  Clearly  the  para¬ 
meter  v  provides  a  measure  of  the  degree  of  randomness  or  "homogeneity"  of 

the  structure.  In  this  sense,  this  model  does  provide  the  pattern  replica- 
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tion  attribute.  For  small  v  the  random  field  (f(x),  xeR  }  appears  as  a 
random  rectangular  mosaic.  As  v  increases,  individual  realizations  rapidly 
approach  a  more  periodic  mosaic  in  appearance.  The  parameters  A,  p  and  v  then 
completely  describe  this  class  of  2-D  randan  fields. 

Although  this  class  of  2-D  random  fields  provides  a  useful  model  of  tex¬ 
ture  in  selected  applications,  the  rectangular  mosaic  exhibited  by  individual 
realizations  is  not  entirely  consistent  with  edge  structure  in  real-world  imagery. 
That  is,  we  would  expect  the  edge  structure  to  exhibit  a  much  more  random  edge 
orientation.  An  alternative  approach  then  is  to  randomly  partition  the  plane 
into  more  complex  geometric  regions.  In  what  follows  we  describe  one  such 
approach  where  the  plane  is  partitioned  into  random  polygonal  regions.  Other 
approaches  are  described  in  [28],  [29]. 

PotygomJL  Pahtlttonl :  Consider  the  partition  of  the  plane  R2  by  a  field  of 
random  sensed  lines.  More  specifically,  an  arbitrary  sensed  line  can  be 
described  in  terms  of  the  3-tuple  (r,0,£  ) .  Here  r  represents  the  perpendic¬ 
ular  or  radial  distance  to  the  line  in  question,  0e[-ir,n]  represents  the 
orientation  of  this  radial  vector,  and  finally  C  is  a  binary  random  variable 
assuming  values  ±1  which  specifies  the  sense  or  direction  imparted  to  this 
line  segment.  The  pertinent  geometry  is  illustrated  in  Fig.  7  for  the  case 
5*1.  By  virtue  of  the  direction  imposed  on  this  line  segment  the  plane  is 
partitioned  into  two  disjoint  regions,  R(right  of  line)  and  L(left  of  line) 
such  that  RUL=R2. 

Now  consider  the  field  of  lines  generated  by  the  sequence  {r^,0^,c^}. 

Here  the  sequence  {r^  represents  the  "event  times"  associated  with  a 


homogeneous  Poisson  process  (N(r),  r>0}  with  intensity  X  events/unit  distance 
evolving  according  to  the  radial  parameter  r.  The  sequence  {0^}  is  i.i.d. 
and  uniform  on  £— 7T ,tt ]  while  is  also  i.i.d.  assuming  the  values  ±1  with 

equal  probability. 

The  field  of  random  lines  so  generated  results  in  a  partition  of  the  plane 
into  disjoint  polygonal  regions.  Gray  levels  are  assigned  as  described  in  [30] 
to  result  in  correlation  coefficient  p  with  gray  levels  in  contiguous  regions. 
Typical  realizations  of  the  resulting  random  field  are  illustrated  in  Fig.  8 

for  selected  values  of  Xg=X/n  and  p.  The  quantity  Xg  represents  the  average 

f 

edge  density  along  any  randomly  chosen  line  segment.  This  random  field  is 
again  described  in  terns  of  the  two  parameters  Xfi,  or  equivalently  X,  and  p. 

This  class  of  2-D  random  fields  cm  be  extended  to  include  more  general  point 
processes  (N(r),  r>0}  controlling  the  radial  evolution;  for  example  *  station¬ 
ary  renewal  processes  with  Gamma  distributed  interarrival  times.  Unfortunate¬ 
ly,  the  analysis  of  the  resulting  processes  becomes  quite  complicated  and  as 
a  result  we  will  not  pursue  this  generalization  here. 

IV.  Second-Order  Properties 

We  turn  now  to  the  second-order  properties  of  the  class  of  2-D  random 
fields  described  in  the  preceding  section.  In  the  interests  of  brevity  the 
treatment  will  be  condensed  and  will  make  extensive  use  of  results  reported 
elsewhere . 

RectonguioA.  PaXtctconA:  As  a  first  step  in  the  development  of  the  covariance 
function,  assume  that  the  randan  orientation  0e[-ir,ir]  has  been  chosen  and  that 
k  transitions  have  occurred  between  the  two  points  x  and  x  +  u  where  we 
assume  for  the  moment  u  >.  0.  It  follows  from  (10)  that 

t  Similarly,  in  the  case  of  rectangular  partitions  it  is  easily  shown  that 
the  average  edge  density  along  any  randomly  chosen  line  segment  is  X  *UX/w. 

We  will  consistently  use  the  subscript  e  to  indicate  edge  density  along  a 
randomly  chosen  line  segment  to  distinguish  from  the  corresponding  unsubscripted 

rate  parameter  X  generating  the  random  partition, 
ft  By  this  we  mean  that  k»kj+k_  where  k. ,  i*l,2,  represents  the  number  of 
transitions  along  each  or  the  orthogonal  axes  which  have  now  been  rotated 
by  0  radians. 
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E{f(x+u)f(x)|0,k}  =  02pk  ;  k=0,l,2 . 

The  conditioning  upon  k  is  easily  removed  according  to 


(14) 


S(f (x+u)f(x) j 0}  =  E  E{f(x+u)f(x) 1 0 ,k)  Pk|0(u)  ,  (15) 


where  pk|0(u)  is  the  probability  of  k  transitions  between  x  and  x+u  given  that 
0  is  acting.  We  exploit  the  stationary  renewal  properties  of  the  point  pro¬ 
cesses  (N^(Z),  i> 0),  1=1,2,  in  writing  this  probability  as  a  function  only  of 
the  displacement  u.  In  particular,  | q (u.)  can  be  evaluated  according  to 


Pkle^  *  ;  k=s0’1—  ’ 


(2),~ 


(16) 


where  q^|g(u^)  is  the  probability  that  (N^(Jl),  Jl>0)  has  undergone  j  transitions 
—  T 

in  the  interval  u^,  i=l,2,  which  depends  upon  u  =(u^,u2)  and  0  according  to 


and 


=  u^cosO+UgSinO, 


u2  =  u^cosS-u^sinO. 


Substituting  (l4)  and  (l6)  into  (15)  we  obtain 


"  .k  r  ,d) 


E{f(jc*»)f(»)|e)«o2  J  p*  I 


(17a) 

(17b) 

(18) 


k=0  4=0 

and  by  simple  rearrangement  of  the  double  summation  in  this  last  expression 
we  find 


E{f(x+u)f(x)|0}=a2  l  l  Pk-J^]|e(\)pJ^^2) 
j=0  k=J 


00  00  I  oi  (19) 

=  a2r  V  „m.(l),~  M.r  V  ~n.(2)  i~ 


Assuming  a  uniform  distribution  for  0,  it  follows  that  the  covariance  function 
becomes 

,  ri t 

(20) 


,  rw 

Rff(x+u,x)  =  ^  I  E{f(x+u)f(x)|0)d0  , 


-IT 
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with  the  integrand  given  by  (19).  While  not  immediately  apparent,  it  is 
easily  shown  that  this  last  expression  depends  only  upon  | |u| |  so  that  the 
resulting  random  field  is  indeed  homogeneous  and  isotropic. 


While  explicit  evaluation  of  (20)  is  in  general,  quite  cumbersome,  it 
can  be  evaluated  in  special  cases.  For  example,  in  the  Poisson  case  V=l, 
it  can  be  shown  [31]  that 

.  2  fir/2 

Rff (I  ImI  I )  *  j  exp{-/2(l-p)A|  |u|  |cos(8-TrA)}d6  »  (21) 

0 

while  the  corresponding  power  spectral  density  computed  according  to  (6) 
becomes 


s  (Si)  = 

fT  fl2+2(l-p)2A2 


i_  r  i 

x2  [n2+(i-p)2A2_ 


(22) 


Typical  covariance  surfaces  together  with  intensity  plots  of  the  corres¬ 
ponding  power  spectral  density  in  the  case  of  periodic  partitions  (i.e.,V“°°) 
are  illustrated  in  Fig.  9-  The  autocorrelation  functions  are  plotted  as  a 
function  of  the  normalized  spatial  variable^  |  |u|  \/l  over  the  range  0<|  |u|  |/£.<3, 
while  the  power  spectral  density  is  plotted  as  a  function  of  the  normalized 
spatial  frequency  variable  fi/2rrA  over  the  range  0<ft/ 2rr A_<5 •  Additional  details 
can  be  found  in  [32].  Explicit  evaluation  of  these  quantitites  for  the  general 
case  of  Gamma  distributed  interarrival  times  is  provided  in  [29]. 

Similarly,  the  conditional  joint  probability  of  f^ffx)  and  f2=f(x+u) 
given  both  the  random  angle  0  and  the  number  of  transitions  k  between  x  and 


x+u  is  easily  shown  to  be  given  by 


tt 


P*fl’f2’-’  £+— l®»k)= 


2rro 

1 

'2TTO 


exp  <- 


exp 


2a2 


|  fJ-2pkfjf2+f2 


20z(l-p2k) 


U(frf2) 


k>0 


k=0 


(23) 


tt 


'Here  l  *  l/A  with  A  the  common  rate  parameter  of  the  two  mutually  independ¬ 
ent  point  processes  which  provide  rectangular  partition  of  the  plane. 

It  is  at  this  point  that  the  Gaussian  assumption  is  crucial. 
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Note  that  this  quantity  is  independent  of  x,  x+u  and  0;  we  will  make  use  of 
this  observation  later. 

The  conditioning  upon  k  in  this  case  is  easily  removed  according  to 


p*fl’f2; 


OO 


-  I 

k=0 


hk(f1,f2)pkj0(u)  , 


(24) 


where  P^jg^u)  has  been  defined  previously  as  the  probability  of  k  transitions 
between  x  and  x+u  given  that  0  is  acting.  We  have  used  h  (f  ,f  )  in  the 
second  expression  of  (24)  in  order  to  emphasize  the  functional  independence 
of  the  spatial  parameters  x  and  x+u  and  the  rotation  angle  0. 

Again  under  the  assumption  of  uniform  distribution  for  0,  the  Joint  p.d.f. 
can  be  evaluated  as 

,  r* 

p*fl*f2’-’  -+-}~  2tt  I  p*fl»f2’-*  2tH!e^d0! 

-IT 


l  hk(fl»f2)pk(ltell).  (25) 

k=o  K  d  * 


where 


Pk( I |u| I ) 


Pk|e(u)d0 


-IT 


(26) 


and  we  have  made  explicit  use  of  the  fact  that  the  integral  on  the  right-hand 
side  of  this  last  expression  depends  only  upon  the  Euclidean  distance  | |u| | . 

It  follows  that  (4)  is  indeed  satisfied  and  hence  the  2-D  random  field  is 
homogeneous  and  isotropic  through  all  second-order  statistics. 

To  complete  the  evaluation  of  the  Joint  p.d.f.  pt^.fgj  |  |u|  | }  it  remains 
to  provide  explicit  evaluation  of  p  ( | |ix| | )  in  (26).  This  has  proven  quite 
cumbersome  in  general,  although  quite  tractable  in  several  important  special 
cases.  For  example,  again  in  the  case  V=1  corresponding  to  Poisson  partitions, 
it  can  be  shown  [29]  that 


i 

i 
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(27) 


’•r'albin11  ^  v 

Pk(  I  Ijil  I )  = - -  I  COS  0  exp{-/2X| |u| |cos0}d0  ;  k=0,l,..  j 

0 

which  does  not  seem  capable  of  further  simplification.  At  any  rate,  this 
expression  is  easily  evaluated  by  numerical  integration.  Substitution  into 
(25)  then  yields  explicit  evaluation  of  ptf^.f^; | |u| | }.  Actually,  for  eval¬ 
uation  and  display  purposes,  it  proves  convenient  to  consider  a  normalized 
version  of  this  joint  p.d. f.  defined  according  to^ 

po{fi.f2il|u||}=c2p{0fi,af2;||u||}  ,  (28) 

which  is  plotted  in  Fig.  10  as  a  function  of  f  ,fg  for  selected  values  of  p 
and  the  normalized  displacement  d'^X^I  |u|  | .  Here  the  point  f^f^O  appears 
in  the  center  and  the  plots  cover  the  range  -3<f£<3,  i=l,2.  Note  the  high 

concentration  of  discrete  probability  mass  along  the  diagonal  f^f^  for  small 
values  of  d'.  This  is  a  direct  result  of  the  high  probability  of  x  and  x+u 
falling  in  the  same  rectangular  regions  and  thus  resulting  in  identical  values 
for  fj=f(x)  and  f2=f(x+u).  This  probability  dimishes  for  increasing  d\  In¬ 
deed, as  indicated  in  Fig.  10,  this  "ridge  line"  along  the  diagonal  has 
virtually  disappeared  for  cP«8.  The  off-diagonal  probability  mass  visible  for 
P  =-0.9  is  a  direct  result  of  the  negative  correlation  while  for  p=0.5»  as 
expected,  there  is  visible  probability  mass  distributed  along  the  main  diagonal. 
For  p*0,  of  course,  this  distribution  is  circularly  symmetric  about  the  origin. 
These  observations  are  more  apparent  in  Fig.  11  which  illustrates  intensity 
plots  of  the  logarithms  of  the  corresponding  p.d.f. 's  in  Fig.  10.  Note,  the 
almost  identical  circularly  symmetric  distributions  which  result  for  large  d' 
independent  of  the  value  of  p.  This  says,  in  particular,  that  if  random 
fields  possessing  different  values  of  p  are  to  be  distinguished  on  the  basis 

The  net  effect  of  this  normalization  is  that  the  f  ,fg  3X33  can  considered 
normalized  to  the  standard  deviation  o. 
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of  second-order  probability  distributions  the  normalized  displacement  d'  must  be 
chosen  judiciously.  Values  of  d'  which  are  either  too  small  or  too  large  provide 
little  discrimination  ability.  Values  of  d*  close  to  unity  appear  to  offer 
the  maximum  discrimination  ability.  We  will  have  more  to  say  on  this  point 
later. 

Polygonal  PaAtttcom :  The  second-order  properties  of  this  2-D  random  field 
have  been  described  in  some  detail  in  [30].  Here  it  is  shown  that,  under  the 
assumption  of  a  Poisson  line  process  generating  the  partitions,  the  auto¬ 
correlation  function  is  given  by 


where  I,  (•)  is  the  modified  Bessel  function  of  the  first  kind  of  order  k. 
k 

Similarly,  the  corresponding  power  spectral  density  is  evaluated  according  to 

Sff(n)  =  "Is7*2-1  j  [  i-2Pcosjtp*  ]  [ ( n/Xg ) 2+( llcos* )2 ] ^  *  (30) 

0 

These  quantities  are  illustrated  in  Fig.  12  for  various  values  of  p.  One  not¬ 
able  characteristic  of  this  random  field  is  that  the  power  spectral  density 
.3/ 

behaves  as  ( Si/X  )  2  for  small  values  of  (n/A  ),  i.e.,  S-Jfi)  has  a  singular- 

e  e  ix 

ity  at  the  origin  except  for  p=-l.  This  high  concentration  of  energy  at  low 
spatial  frequencies  is  a  direct  result  of  the  construction  procedure  which 
allows  relatively  large  correlations  between  gray  levels  in  regions  relative¬ 
ly  far  apart.  We  feel  that  this  characteristic  is  typical  of  selected  texture 
processes  and  as  a  result  it  was  purposely  built  into  the  construction  proced¬ 
ure. 

Finally,  following  the  procedure  described  previously  in  the  case  of 
rectangular  partitions,  the  Joint  p.d.f.  is  easily  shown  to  be  given  by  (25) 
with  the  sum  extended  over  both  positive  and  negative  values  of  k  and 
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ik|  *  U_  I  lil  |/2]2*' 

Pk(||«l|)  s  txJ  |u|  |/2]|k|exp{-Aet  |u|  |}  l  {jtV|k]TiTi  ;  k=°.±1’±2»-- (30) 

it— 0 

In  Fig.  13  we  provide  intensity  plots  of  the  logarithm  of  | ju| | )  as  a 

function  of  f^  and  for  selected  values  of  p  and  d'3Xel IhJ  I •  An  inter¬ 
esting  observation  to  be  drawn  here  is  the  persistence  of  the  diagonal  "ridge 
line"  with  increasing  values  of  d'.  This  is,  of  course,  a  direct  result  of 
the  construction  procedure  which  allows  return  to  the  same  gray  level  at 
distant  spatial  locations  with  relatively  high  probability. 


Coirment:  The  comments  made  previously  in  relationship  to  Fig.'s  1  and  2  can 
now  be  explained  in  terms  of  the  preceding  properties  of  the  2-D  random  field 
models  of  texture.  For  example,  in  the  case  of  rectangular  Poisson  partitions 
the  autocorrelation  function  and/or  power  spectral  density*  computed  according 
to  (21)  and  (22)  respectively* depend  only  upon  the  product  (l-p)x  .  The 
sample  fields  in  Fig.  1  have  been  chosen  to  maintain  a  constant  value  for 
this  product  and  hence  possess  identical  second-moment  properties.  Similar¬ 
ly,  for  this  random  field  the  average  edge  density  along  any  randomly  chosen 
line  segment  can  be  shown  to  be  given  by  A  =**A/ir.  The  sample  random  fields 
in  Fig.  2  all  possess  identical  values  for  Xg  and  hence  cannot  be  distinguish¬ 
ed  on  the  basis  of  either  edge  density  or  average  gray-level  run  lengths 
alone.  Finally,  as  we  have  observed  in  conjunction  with  Fig.’s  10  and  11, 
it  is  possible  for  two  visually  distinct  random  fields  to  possess  the  same 
joint  p.d. f.  for  some,  but  clearly  not  all,  displacement  distances  d'.  Again 
this  points  out  a  potential  problem  associated  with  discrimination  approaches 
based  upon  spatial  gray-level  co-occurrence  probabilities. 
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V .  Log-Likelihood  Texture  Discriminator: 

N 

We  suppose  that  each  point  of  an  image  array  {f.  }  is  known  to  be 

i » J  i » J=i 

associated  with  one  of  a  finite  number  K  of  texture  classes  or  hypotheses 

N 

labeled  H. ,  i=0,l, — ,K-1,  respectively.  Furthermore,  the  array  {f.  }. 

i  i ,  J  i ,  1 

will  be  assumed  to  have  been  obtained  by  discrete  homogeneous  sampling  of 

2 

corresponding  sample  functions  (f(x),  x£^  }  of  2-D  random  fields  as  described 
in  preceding  sections.  Typical  realizations  are  illustrated  in  Fig.  lL  in 
che  case  of  the  rectangular  partition  process.  Here  various  texture  samples 
are  provided  in  different  corners  of  an  image  which  are  labeled  as  points  on 
a  compass,  i.e.,  NW  indicates  the  northwest  corner,  E  indicates  the  east,  etc. 
In  the  top  row  of  Fig.  1^,  the  parameter  v  is  held  fixed  for  each  image  while 
p  is  varied.  Similarly,  in  the  bottom  row  p  is  held  fixed  for  each  image 
while  v  is  varied.  In  all  cases  the  edge  density  is  held  fixed  at  A=0.50. 
Additional  examples  have  been  provided  previously  in  Fig.'s  1  and  2.  Our 
interest  then  will  be  in  discriminating  between  the  various  texture  regions 
and  accurately  detecting  the  boundaries  between  these  regions. 

An  optimum  texture  discriminator,  which  minimizes  the  classification  error 
should  clearly  be  based  on  a  threshold  test  on  the  likelihood  functional  or 
some  monotonic  function  of  it  [33],  [3^].  More  specifically,  suppose  a  window 
of  size  (2M+l)x(2M+l)  is  constructed  about  each  pixel  position  (i,j).  The 
observations  at  pixel  positions  within  the  window  centered  at  position  (i»j) 
will  be  denoted 

Fi,j  =  {fk,A;  5  i,j=l,2,...,N  ,  (32) 

where  we  assume  M<<N  and  neglect  the  boundary  effects.  The  likelihood 
functional  with  the  window  centered  on  pixel  position  (i,j)  and  assuming  the 


19 


i 


k'th  hypothesis  acting  is  then  defined  according  to 


a  P<F,  A1 

MF<  ,isritT  1  k»o,i,...,K-i  . 
*  1>J  po{Fi,j} 


(33) 


where  Pq(*)  is  an  appropriately  defined  p.d.f.  independent  of  which  hypothesis 

is  acting  and  serving  merely  to  provide  a  convenient  normalization.  As  the 

N 

window  is  scanned  accross  the  image  array  {f.  }.  ,  each  pixel  position 

i  >  J  i J 

is  assigned  the  value  k(i,j)  corresponding  to  the  index  which  maximizes  the 
class-conditional  likelihood  functional  in  (33).  That  is,  k(i,j)  =  k^  if 


\  =  max  \{Fi 

0  0<k<K-l 

In  this  way  each  point  will  be  assigned  to  one  of  a  predetermined  collection 

of  texture  classes  on  the  basis  of  local  gray  level  values  in  a  way  that 

minimizes  the  classification  error.  Boundaries  can  be  determined  by  edge 

N 

detection  of  the  resulting  array  {k(i,j)}.  .  .  . 

l  *  j— J- 

There  are  several  difficulties  with  the  approach  described  above,  first, 
in  order  to  implement  the  likelihood  functional  in  (33)  we  require  explicit 
knowledge  of  the  joint  p.d.f.  of  order  (2M+l)x(2M+l)  for  the  underlying  2-D 
random  field  (f(x),  xeR  }.  It  is  very  rare  that  this  information  would  be 
available.  For  the  2-D  random  fields  described  in  the  preceding  sections 
we  have  been  able  to  determine  the  joint  p.d.f.  ptf^.f^;!  |u|  |}  for  two  points 
separated  by  the  distance  | |u| | ;  higher  order  p.d.f.'s  are  somewhat  intractable. 
Furthermore,  even  if  these  high-order  p.d.f.'s  were  available  any  texture  dis¬ 
crimination  approach  based  upon  the  likelihood  functional  would  be  so  finely 
tuned  to  the  modeling  assumptions  that  it  may  be  of  questionable  utility  in 
a  practical  application.  We  might  prefer  a  suboptimum  al^ though,  hopefully, 
more  robust  approach. 
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A  reasonable  suboptimum  approach  would  be  to  perform  a  data  reduction 

on  the  sequence  (F.  .}  of  observations  to  result  in  measurements  which  are 

i  .J 

more  tractable  um’  at  the  same  time  relatively  robust.  The  reduced  data 

can  then  be  used  in  a  subsequent  log-likelihood  discriminator.  There  are 

many  candidates  for  the  data  reduction  operation.  We  have  found  it  useful 

to  make  use  of  histograms  of  gray-level  co-occurrences.  More  specifically, 

suppose  that  the  original  image  array  {f.  }.  is  quantized  to  Q  levels 

labeled  0,1,..., Q-l.  Haralick  and  his  co-workers  [ 1 ]— [ 2 ]  have  introduced 

the  concept  of  the  gray-level  co-occurrence  matrix.  This  is  a  QxQ  matrix 

with  (m,n)  element  p.  {m,n;  d,0}  defined  as  the  number  of  times  the  gray 

1  >  J 

levels  m  and  n  occur  separated  by  d  pixels  at  an  angle  0  within  a  window  of 

size  (2M+l)x(2M+l)  centered  at  pixel  position  (i,j).  Furthermore,  since 

the  texture  models  under  consideration  are  invariant  under  rotations  we 

utilize  a  version  of  the  co-occurrence  matrix  which  is  an  average  over  all 

0  and  denote  the  result  by  P.  (d).  The  proposed  texture  discriminator  is 

1  >  J 

then  implemented  as  a  log-likelihood  test  on  the  sequence  {p.  (d)}  for  a 

1  »  J 

fixed  value  of  d.  A  block  diagram  of  the  resulting  log-likelihood  dis¬ 
criminator  is  illustrated  in  Fig.  15.  The  pixel  position  (i,j)  is  assigned 
the  value  k(i,j)=kQ  if 


Lfc  (P,  *(<*)}  =  max  L{p  ,(d)}  , 

0  0<k<K-l 


(35) 


M)\tV 

P  TdH~  ;  k=0,l,...,K-l  ,  (36) 

i » j 

represents  the  class -conditional  log-likelihood  functional. 

Haralick  has  advocated  use  of  the  gray-level  co-occurrence  matrix 
differently  than  that  described  here.  More  specifically,  he  has  proposed 


p(P 


Po' 
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lU  functionals  or  discriminants  defined  on  the  gray-level  co-occurrence 
matrix  which  are  intended  to  discriminate  among  different  texture  regions. 
While  these  quantities  have  been  shown  to  be  of  value  in  selected  applica¬ 
tions  they  are  based  upon  heuristics  for  the  most  part.  Here  the  use  of 
P.  ,(d)  in  a  maximum  likelihood  detector  is  more  solidly  based  on  statistical 

decision  theory  concepts.  The  primary  motivation  for  employing  P.  (d)  to 

1  >  J 

effect  the  data  reduction  in  Fig.  15  is  the  analytical  simplicity  which 

results.  The  conditional  p.d.f.'s  p{P.  (d)|H  }  ,  k=0,l,. . . ,K-1  are  easily 

evaluated  in  terms  of  the  joint  p.d.f.  p{f1>f^; | |u| | }as  demonstrated  below. 

Use  of  second-order  p.d.f.'s  provides  a  logical  extension  of  techniques 

based  upon  either  correlation  properties  or  edge  density.  This  choice  has 

also  been  influenced  by  the  general  inability  of  humans  to  discriminate 

textures  which  differ  only  in  higher  order  statistics. 

In  the  implementation  of  the  log-likelihood  texture  discriminator  a 

critical  simplifying  assumption  is  made  concerning  the  construction  of  the 

P  (d)  matrix.  Specifically,  we  assume  that  each  pair  of  points  counted  in 
*  »  J 

producing  P  (d)  has  its  intensities  assigned  independently  of  every  other 
*  »  J 

pair  of  points.  Clearly  this  is  not  generally  the  case  for  the  2-D  random 
fields  described  in  preceding  sections.  Nevertheless,  justification  to 
some  extent  follows  from  the  assumed  ergodic  properties  of  the  underlying 
random  fields.  In  particular,  we  assume  that  spatial  averages  over  the 
observation  window  approach  ensemble  averages  asymptotically  as  the  size  of 
the  window  increases.  As  a  result,  the  empirical  distributions  constitut¬ 
ing  elements  of  p.  ,(d)  approximate  the  true  probability  distributions  of 
1  *  J 

corresponding  spatial  gray-level  co-occurrence  events.  Hence,  provided  the 
observation  window  is  large  enough,  the  matrix  P.  ,(d)  is  indistinguishable 
from  that  which  would  have  been  obtained  from  a  set  of  independent  samples 
of  the  same  process. 
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The  description  of  the  log-likelihood  texture  discriminator  illustrated  in  Fig.  15 
is  then  complete  once  explicit  evaluation  of  the  class-conditional  log-likelihood 
functionals  L  {P.  (a)}  is  provided.  Under  the  independence  assumption  de- 

scribed  above,  the  evaluation  of  p{P  (d)|K  }  ,  k=0,l, . . . ,K-1,  is  described 

A  »  J  & 

in  terms  of  the  multinomial  distribution  [15], 

Q-l 


/  yrl  \ 

— J--H 


Q-l 

TT 

m,n=0 


(P,  ,(m,n;d)l)m’nS° 


i.J 


p  , (m,n;d) 
1  1 


Vn<d;V 


(37) 


where  n(d;H^)  is  the  probability  of  observing  gray  levels  m  and  n  at  a 

particular  pair  of  points  separated  by  distance  d  under  hypothesis  H^,  k=0,l, 

. . . ,K-1.  For  convenience  we  take  the  normalization  functional^ p  { P.  (d)} 

0  1 » J 

/  Q-l  \ 

1  l  p  (m,n;d))l 


as 


FT  (Pj^fm.nid)!) 


(38) 


m,n=0 

It  follows  from  (36)  that  the  class- conditional  log-likelihood  functionals  are  given 
by  Q-l  Q-l 

VPi  (d)}  =  l  l  p  (m,n;d)JUm,n;d)  ;  k=0,l, . . .  ,K-1  ,  (39) 

*  1,J  m=0  n=0  k 

where 

«.k(m,n;d)  =  In  ^  n(d;Hk)  •  (UO) 

This  latter  quantity  is  easily  computed  from  the  joint  p.d.f.  of  the  underly¬ 
ing  2-D  random  field.  Suppose  that  the  {fk(x),  xeR^)  represents  the  random 
field  associated  with  hypothesis  class  k=0,l, . . . ,K-1,  and  possessing 


joint  p.d.f.  Pk<f1*f2;l  lil  !>•  711611 

f^m+1  rEn+l 

Vn(dlV  *  J  J  >>k<fl-f2‘d4I-)dfldf2  • 


Hi) 


m 


Strictly  speaking  this  is  not  a  p.d.f.  but  does  provide  a  convenient  normal¬ 
ization  of  the  class-conditional  log-likelihood  functionals. 
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where  E^,  £=0,1,. . . ,Q-1  represents  the  lower  boundaries  of  the  quantization 
bins  associated  with  level  £,  and  AL  represents  the  spatial  sampling  interval. 

The  likelihood  functional  defined  by  (39)  has  a  simple  interpretation  as 
an  inner  product  operation.  More  specifically,  let  L^(d)  represent  the  QxQ 
array  with  (m,n)  element  £^(m,n;d).  Then  (39)  becomes 


Lk{Pi,j(d)}  =  <Pi,j(d)»  Lk.(d)>  '•  k=0,l,...,K-l 

e  use  of  the  i 

<A,  8>=  tr  AB1 


where  we  have  made  use  of  the  inner  product 


(1*2) 

(1*3) 


defined  on  the  space  of  QxQ  matrices.  The  precomputed  elements  of  Lfc(d)  can 

be  stored  in  a  table  and  accessed  as  required.  Computation  of  the  inner 

product  in  (1*2)  is  relatively  simple.  The  major  computational  burden  is 

associated  with  computation  of  the  P.  (d)  matrix.  In  the  next  section  we 

£ » J 

describe  a  simple  digital  implementation  of  the  log-likelihood  discriminator 

which  does  not  require  explicit  computation  of  P.  ,(d).  This  approach  is 

1  ♦  J 

based  upon  2-D  Weiner  filtering  concepts.  The  result  is  a  computationally 
efficient  implementation  as  a  2-D  recursive  or  infinite  impulse  response 
(IID)  digital  filter.  Furthermore ,  the  recursive  implementation  leads  to  a 
simple  method  of  avoiding  the  block  structure  associated  with  the  windowing 
operation. 


VI.  Digital  Filter  Implementation  of  Log- Likelihood  Discriminator: 

The  log-likelihood  functional  described  by  (39)  and  the  sequel  requires 
a  summation  over  gray  levels  or  intensities.  This  can  be  replaced  by  a  spatial 
summation  and  leads  to  a  simple  implementation  as  a  2-D  digital  filtering 
operation.  In  particular,  recall  that  p.  , (m,n;d)  is  simply  the  number  of 

*  »  «J 

times  that  gray  levels  m  and  n  occur  at  a  pair  of  points  separated  by  d 
pixels  within  the  window  W,  . ,  centered  at  pixel  position  (i,j).  It  follows 


2k 


that  the  class- conditional  log-likelihood  functional  can  be  rewritten 
as 


VFi  wa>  =  l  l 

*  (r,s )tW. 

*  J 


l  I  Mf  c+vid) 

(u,v)e5(r,s;d)  AWi^  r’s  r  u’s  v 

»k=0,l, . . . ,K-1, 


(U4) 


where  S(r,s;d)  is  the  set  of  all  points  which  are  distance  d  from  pixel 

position  (r,s).  The  outer  summation  in  (UU)  is  over  the  pixel  positions  (r,s) 

within  the  window  W.  centered  at  (i,j)  while  the  inner  summation  is  over  the 
1  *  J 

set  S(r,s;d)fiW.  .  including  only  those  points  (u,v)  within  the  window  which 

1  9  J 

simultaneously  stand  in  the  specified  spatial  relationship  to  pixel  position 
(r,s).  Note  that  this  latter  summation  can  be  obtained  by  searching  over  a 
circular  neighborhood  of  radius  of  at  most  d  units.  We  assume  that  appro¬ 
priate  bookkeeping  has  been  employed  to  avoid  double  counting  of  pairs  of 
points  in  the  specified  spatial  relationship.  This  is  indicated  by  the  prime 
added  to  the  inner  spatial  summation  in  (kk) . 

At  this  point,  it  is  convenient  to  define  the  quantity 


g^(r,s;d)  =  l  •  l'  s’fr+u  s+v;d)  i  ^=0,1, . . . ,K-1. 

^  (u,v)eS(r,s;d) fi  W  .  k  r’s  r+u’s+v 

1  »  J 


(^5) 


We  are  justified  in  expressing  gv(r,s;d)  as  a  function  only  of  spatial 


coordinates  (r,s)  since  once  the  sampled  field  {f.  has  been  observed 

1  >  J  1  »  J  J- 

this  quantity  is  readily  evaluated  through  (1*0)  and  the  sequel.  This  assumes, 

of  course,  that  boundary  effects  along  the  periphery  of  the  window  W.  can 

be  neglected;  otherwise  g^(r,s;d)  would  likewise  depend  upon  the  spatial 

index  ( i , j ) .  These  boundary  effects  will  be  assumed  negligible  in  what 

follows.  This  is  a  reasonable  assumption  if  d<<M  (recall  the  window  w.  ,  is 

1  *  J 

of  size  ( 2M+l)x(2M+l) )  which  will  generally  be  the  case.  It  follows  that 
(M)  becomes 
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Lk{i,J;d>  -  l  l  g^r.sjd)  ;  k*0,l,. • .  .K-l  ,  (U6) 

where  we  have  found  it  convenient  to  suppress  the  functional  dependence  upon 
the  observations  and  simply  write  L^Ci.Jid}  for  jjd). 

Consider  now  how  the  window  size  might  be  expected  to  influence  perform¬ 
ance  of  the  log-likelihood  discriminator.  Clearly  a  large  window  size  is 
desirable  if  the  window  covers  an  area  of  homogeneous  texture,  as  this  re¬ 
duces  the  probability  of  classification  error.  As  the  window  size  is  increas¬ 
ed,  however,  it  becomes  more  likely  that  the  window  will  contain  two  or  more 
regions  of  different  texture  classes  thereby  increasing  the  classification 
error.  To  arrive  at  a  reasonable  compromise  between  these  conflicting  factors 
the  window  has  been  allowed  to  be  of  infinite  extent^but  the  sum  in  (1*6)  is 
replaced  by  a  weighted  simulation.  The  weights  can  then  be  chosen  to  provide 
diminishingly  less  weight  to  points  the  further  they  are  from  pixel  position 
(i,j).  More  specifically,  we  generalize  (U6)  to  the  form 

h(i-r,J-s)gk(r,s;d)  ;  k=0,l,. . .  ,K-1  ,  (1*7) 

r  s 

with  (h(r,s)}  the  weighting  function.  Note  that  (1*7)  reduces  to  (1*6)  under 

the  assumption  . 

Jl  ;  | r [ , | s  j <M 

h(r,s )  *  < 

lo  ;  elsewhere  ,  (1*8) 

Furthermore,  observe  from  (1*7)  that  I^ti.jjd}  is  simply  the  output  of  a  2-D 
digital  filter  with  the  sequence  {g  (r,s;d))  as  input.  The  weighting  sequence 

K 

(h(r,s)}  is  the  2-D  impulse  response  cr  point  spread  function  of  this  filter. 


Note  that  this  justifies  our  assumptions  that  boundary  conditions  are 
negligible  in  computing  gk(r,s;d)  from  (1*5). 
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There  are  many  heuristic  arguments  that  can  be  developed  for  choosing 
an  impulse  response  sequence  {h(r,s)}.  We  have  found  it  useful  to  formulate 
this  choice  as  a  2-D  Wiener  filtering  problem.  More  specifically,  consider 
the  filter  input  sequence  defined  by  (45)  as  a  function  of  a  continuous 

p 

spatial  variable  x.  The  quantity  {g^(x;d),  xeR  }  is  then  a  2-D  random 
field  for  each  k=0,l, . . . ,K-1,  and  parameterized  by  the  distance  d.  Although 
the  exact  nature  of  this  field  is  rather  difficult  to  describe  precisely, 
we  will  make  some  rather  crude  modeling  assumptions  which  have  led  to  some 
useful  results.  In  particular,  we  assume  that  g^(x;d)  can  be  represented 
as  the  sum  of  three  separate  components  according  to 

g^xid)  *  t(x)+  i(x)  +n(x)  ;  xeR2  .  (49) 

Here  t(x)  represents  a  mean- value  or  signal  component  indicative  of  the  true 
texture  class.  This  component  is  constant  over  homogeneous  texture  regions 
and  exhibits  Jumps  or  discontinuities  at  the  boundaries  between  different 
textured  regions.  The  component  i(x)  represents  an  interference  component 
the  nature  of  which  is  similar  to  the  texture  within  a  region.  This  compon¬ 
ent  reflects  the  pixel-to-pixel  variations  in  g^xjd)  due  to  residual  texture 
components.  Finally,  n(x)  represents  an  unavoidable  noise  component  repre¬ 
senting  background  noise,  spurious  image  detail,  quantization  noise,  etc. 

In  what  follows  the  signal  component  t(x)  will  be  modeled  in  terms  of  a 
polygonal  partition  process  as  described  in  Section  III.  The  edge  density^  *te> 
or  equivalently  X^,  will  be  chosen  on  the  basis  of  an  assumed  density  for 
texture  boundaries  while  the  correlation  p  =0  will  be  chosen  to  reflect  com- 
plete  independence  of  texture  in  contiguous  regions.  Note  that  this  latter 
assumption  implies  that  the  distance  d  has  been  chosen  appropriately  to 
maximize  the  discrimination  ability;  otherwise  the  signal  process  t(x) 

t  Again,  by  our  subscripting  convention  X  represents  the  edge  density  of  the  field 
t(x)  along  a  randomly  chosen  line  segment. 


would  be  highly  correlated  in  contiguous  textured  regions.  The  interference 
process  i(x:),  on.  the  other  hand,  will  be  modeled  as  a  texture  process,  either 
rectangular  or  polygonal,  possessing  much  higher  edge  density  and  with  para¬ 
meters  chosen  to  match  the  coarsest  texture  expected  in  the  input  image.  In 
a  sense  this  represents  a  worst  case  choice  as  the  coarsest  (i.e.,  lowest 
edge  density)  texture  presents  the  most  difficulty  in  separation  from  t(x) 
by  linear  filtering.  Finally,  the  noise  field  n(x)  is  assumed  a  white  noise 
field  with  power  spectral  density  S^CQ)  =a^. 

It  should  be  noted  that  the  model  specified  by  ( 1+9 )  is  independent  of 
both  k  and  d,  whereas  in  reality  we  would  not  expect  this  to  be  the  case. 
Nevertheless  this  model  has  led  to  some  useful  suad  interesting  results  in 
selected  computer  experiments,  some  of  which  are  described  in  the  next 
section.  The  main  Justification  for  these  modeling  assumptions  is  based  upon 
rather  extensive  empirical  observations  on  typical  realizations  of  {gk(x;d),xeR  }. 
Although  the  modeling  is  somewhat  crude,  it  does  provide  some  consideration 
of  the  relationship  between  the  size  of  regions  of  constant  texture  and  the 
coarseness  of  this  texture  in  the  design  of  the  filtering  operation.  Spec¬ 
ifically,  there  is  an  inherent  tradeoff  between  the  degree  of  smoothing  of 
point-to-point  variations  of  gjJxjd)  and  the  ability  of  the  generalized  log- 
likelihood  discriminator  to  distinguish  small  regions  of  homogeneous  texture 


from  contiguous  regions.  The  Wiener  filtering  problem  then  is  to  determine 
the  linear  least  mean-square  estimate  of  t(x)  from  the  noisy  observations 
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where  and  S^(fi)  are  the  power  spectral  densities  of  the  true  texture 

component  t(x)  and  the  point-to-point  interference  component  i(x)  respectively. 
As  previously  demonstrated,  under  our  modeling  assumptions,  these  quantities 
depend  only  upon  the  radial  frequency  Q  and  hence  the  Wiener  filter  possesses 
this  symmetry  property.  The  system  transfer  function  in  (56)  can  be  rewritten 
as 


H0(n)  = 


1 

s...  (a) 


1+ 


ii 


tt 


,ys) 

m  stt 


(n) 


(52) 


which  suggests  defining  the  parameters  y^  and  y^  as  the  ratio  of  interference 
and  noise  powers  to  the  power  in  the  texture  process  respectively.  These 
parameters  provide  a  measure  of  the  degree  of  degradation  by  interference  from 
residual  texture  components  and  salt-and-pepper  noise.  More  specifically, 

(53a) 
(53b) 

so  that  the  Wiener  filter  transfer  function  HQ(ft)  is  completely  determined  in 
terms  of  the  quantities  y^ >Vn»^ and^  p^.  In  any  particular  application 
these  can  be  estimated  empirically  or  on  the  basis  of  a  priori  knowledge  con¬ 
cerning  the  texture  classes  to  be  discriminated. 

The  preceeding  description  of  the  optimum  linear  filtering  operation  is 
for  a  continuous  spatial  domain;  we  have  access  only  to  sampled  data  so  that 
a  digital  implementation  is  required.  In  particular,  we  seek  a  2-D  digital 
filter  with  system  transfer  function  H^z^z,,)  whose  frequency  response++ 
approximates  HQ(fl).  While  there  are  many  digital  implementations  available, 

*  Actually  Pt=0  in  the  work  described  here. 

tt  The  frequency  response  of  the  2-D  digital  filter  is  simply  H  (z^,z„)  evaluated 
for  z  =eJ  l,  i=l,2.  We  assume  the  spatial  frequency  variable  Is  measured 
in  units  of  radius  per  sample  distance. 


we  have 


while 


Yi  ’'IK 


Y  =  a2/a2 

yn  un/  t 
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our  approach  has  been  to  make  use  of  a  2-D  infinite  impulse  response  (IIR) 
digital  filter  whose  point  spread  function  exhibits  four-quadrant  symmetry. 

The  choice  of  an  IIR  filter  was  dictated  by  the  computational  economies  which 
result  from  the  inherent  recursive  nature  of  the  computations  involved.  Simil¬ 
arly,  the  restriction  to  four-quadrant  symmetry  for  (h(r,s)}  follows  from  the 
requirement  that  the  log-likelihood  discriminator  not  exhibit  any  directional 
sensitivities. 

The  IIR  filters  will  be  assumed  to  possess  rational  system  transfer  functions 


of  the  form 


*b  Nt 


VVV  =  l  I  bii 

i=0  j=0 


'  M  N 

a  a 

1+  £  E  aij  ZI  Z2J 

i=0  j=0  J 

i-J*0 


In  particular,  the  output  sequence  {y.  , }  in  response  to  {x.  .}  as  input  can 

i ,  J  i » J 

be  obtained  recursively  according  to 


M  N  M.  N 

a  a  d  b 

/_  =  -  l  I  a  y  .  „  .+  I  l  b  x  ,  ;  m,n>0  •  (55) 

i=0  J*0  ^  i*n—  J  q  j=Q  m— i  ,n— j 

i=J*) 


We  will  assume  that  the  geometry  is  such  that  this  corresponds  to  a  stable 
filter  recursing  from  the  upper  left-hand  corner.  Observe  that  the  result¬ 
ing  filter  will  have  nonzero  impulse  response  only  in  the  lower  right  quadrant. 

As  mentioned  previously,  it  is  highly  desirable  that  the  digital  filter 
implementation  of  the  log-likelihood  discriminator  exhibit  four-quadrant 
symmetry  in  its  impulse  response  or  point  spread  function.  One  method  for 
achieving  a  point  spread  function  with  this  inherent  symmetry  is  to  allow 
repeated  application  of  the  same  filter  recursing  from  each  of  the  four  corners. 
If  {h  (i,j)}  represents  the  point  spread  function  associated  with  a  single 


application  of  the  filter  specified  by  (5*0  then  the  composite  filter  possesses 


point  spread  function  as  indicated  in  Fig.  16.  The  corresponding  system  trans¬ 
fer  function  of  the  composite  filter  is  then 


H(z1>z2)  =  z^  ^Zg”  ^H^z^Zg)  , 


(56) 


where 


H1^Z1’Z2^H0^Z1,Z2^+Z1H0^Z1~1’Z2^+Z2H0^Z1,Z2~1^+Z1Z2H0^Zi”1,Z2~1^' 

It  follows  that  the  frequency  response  of  H^.z^)  is  identical  to  that  of 

H^( z^,Zg)  up  to  an  unimportant  linear  phase  term.  In  choosing  YLAz  to 

provide  an  approximation  to  H An)  for  z  =e^Wi,  i=l,2,  we  have  restricted 

0  1 

attention  to  the  case  where  Hq(z^,z2)  is  a  simple  first-order  section  of  the 
form 

-1, 


1'"bio(zi1't'z2i)'t‘1:>nzilz21 
1+a10(*I1+*21)',,alIZilB21 


H0(Z1,Z2)=A 


(58) 


This  choice  insures  zero  frequency  response  at  the  origin  and  symmetry  of  the 
corresponding  point  spread  function  about  a  line  at  U50  to  the  axes,  i.e., 
hQ(i,j )=hQ( j ,i) .  Similar  properties  extend,  of  course,  to  the  composite  filter 
represented  by  H(z^tz^).  A  computer  program  has  been  written  for  determining 
the  four  coefficients  a10,an’  ^lO’^ll  ga*n  A  according  to  an  itera¬ 

tive  gradient  procedure  to  result  in  a  frequency  response  for  H^(z^,z2)  which 
provides  a  least  mean-square  approximation  to  the  desired  response  HQ(S2).  The 
details  of  this  program  are  described  in  [29]. 

In  Table  I  we  summarize  the  results  of  this  iterative  digital  filter  design 

approach  for  selected  values  of  y.,p.,A.,p  ,  and  A  all  with  y  =  -lOdB.  Here 

x  1  1  v  w  n 

we  have  found  it  convenient  to  classify  the  interference  characteristics  as 
weak,  moderate  or  strong  depending  upon  the  value  of  y^.  The  parameters  chosen 
here  are  particularly  relevant  to  some  experimental  results  to  be  described  in 
the  next  section. 
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Interference 

Characteristic 

Xi 

pi 

Xt 

pt 

A 

— 

a,  „ 

10 

all 

bio 

bll 

Weak 

0.16 

0.0 

0.0125 

0.0 

0.8131 

-0.97*+3 

0.952U 

0.2915 

-0.9098 

Moderate 

0.16 

0.0 

0.0125 

0.0 

0.2032 

-0.9^93 

0.9036 

-0.0375 

-0.1761 

Strong 

0.16 

0.5 

0.0125 

0.0 

0.U220 

-0.9605 

0.9220 

-0.8796 

0.8U32 

Table  1 

Typical  Filter  Parameters  for  Yn=  -lOdB 

For  the  three  cases  described  in  Table  1,  the  corresponding  power  spectral 

densities  S^(ft)  and  S..(n)  +  S  (Q),  together  with  the  resulting  Wiener  filter 
tt  li  nn 

response  HQ(fi),are  plotted  in  Fig.  17  as  a  function  of  the  radial  frequency 
variable  n.  Observe  the  lowpass  behavior  in  all  cases  with  the  selectivity 
increasing  with  increasing  levels  of  interference.  Finally,  in  Fig.  18  we 
illustrate  3-D  plots  of  the  desired  Wiener  filters  and  the  resulting  digital 
approximations.  In  all  cases,  the  closest  corner  represents  the  point^  [-ir,ir] 
while  the  farthest  corner  represents  the  point  The  left-hand  column 

shows  the  desired  responses  while  the  right-hand  column  illustrates  the  re¬ 
sponses  exhibited  by  the  digital  approximations.  Note  that  in  all  cases  the 
lowpass  nature  of  the  optimum  filter  has  been  preserved  and  a  fair  degree  of 
symmetry  has  been  retained. 


t  Frequency  axes  for  w, ,  1*1,2,  have  been  normalized  to  the  range  [-n,ir]  with 
the  endpoints  corresponding  to  ±  the  folding  frequency,  i.e.,  half  the 
sampling  rate. 
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VII. 


Experimental  Results : 

Typical  performance  results  for  selected  texture  discriminators  are 
provided  in  Fig.  19.  Here  Fig.  19a  illustrates  the  original  image  which 
consists  of  realizations  of  three  distinct  rectangular  partition  processes 
for  various  parameter  choices.  The  NW  and  NE  corners  have  A=0.l6,  p=0.0, 
and  A=0.32,  p=0.5  respectively.  These  values  were  carefully  chosen,  accord¬ 
ing  to  previous  comments,  to  result  in  identical  second-moment  properties. 

As  a  result  these  two  fields  cannot  he  discriminated  on  the  basis  of  autocorrel¬ 
ation  functions  and/or  power  spectral  densities  alone.  The  field  in  the  S 
corner  has  A=0.32  while  p=0.0.  Since  it  possesses  the  same  edge  density  as 
the  field  in  the  NE  corner,  these  two  textures  cannot  be  discriminated  on  the 
basis  of  edge  density  alone.  Although  somewhat  copt rived,  we  feel  that  this 
problem  provides  a  real  challange  to  texture  discrimination  algorithms. 

Three  log-likelihood  discriminators  were  designed  using  the  filter  para¬ 
meters  corresponding  to  the  three  entries  in  Table  1.  The  choice  A .=0.16 
corresponds  to  the  coarsest  texture  in  Fig.  19a  while  A  was  chosen  to  approximate 
the  density  of  texture  boundaries.  Recall  that  the  edge  density  in  this  case 
is  given  by  A.  =A. /ir  edges  per  pixel  so  that  for  this  choice  there  would  be  on 
average  approximately  one  texture  transition  along  the  boundary  of  256x256  image. 
This  approximates  the  situation  illustrated  in  Fig.  19a. 

In  Fig.  19b  we  illustrate  the  performance* of  the  log-likelihood  texture 
discriminator  for  the  case  of  moderate  interference,  i.e.,  the  middle  entry 
in  Table  1.  The  value  of  d  used  here  was  d=3.  Ideally,  this  quantity  should 
be  chosen  such  that  &*l/Ae  (recall  Ae=l*A/it  for  the  rectangular  partition 
process)  to  provide  maximum  discrimination  ability.  In  situations,  such  as 
that  presented  in  Fig.  19a,  where  there  are  more  than  one  value  of  edge  density 

t  In  all  experimental  results  reported  in  this  section  the  images  were 
uniformly  quantized  to  Q=64  levels. 


associated  with  the  texture  samples  to  be  discriminated,  a  reasonable  choice 

is  to  choose  d  *lAe  where  Ag  represents  the  average  edge  density  over  all 

texture  classes.  For  the  present  case  we  have  A  =  .306  and  hence  the  choice 

e 

d=3  provides  the  desired  approximation.  Studies  have  indicated  that  perform¬ 
ance  is  not  a  sensitive  function  of  d  in  the  range  2<d<3, although  the  choice 
d=3  appeared  to  be  about  optimum.  As  indicated  by  Fig.  19b,  the  log-likeli¬ 
hood  discriminator  does  an  excellent  job  of  discriminating  the  three  texture 
regions  except  in  the  vicinity  of  either  texture  or  image  boundaries.  This 
performance,  however,  can  be  improved  as  we  demonstrate  subsequently. 

Included  in  Fig.  19  for  comparison  purposes  we  have  indicated  the  perform¬ 
ance  of  alternative  more  conventional  texture  discrimination  schemes.  In  Fig. 

19c  we  demonstrate  the  performance  of  a  conventional  correlation  discriminator. 

This  algorithm  implements  a  threshold  test  on  a  least-squares  estimate  of  the  correl¬ 
ation  of  pixels  separated  by  distance  d.  The  optimum  threshold  has  been  chosen 
empirically  on  the  basis  of  histogram  techniques.  While  this  approach  is  useful 
in  discriminating  the  texture  in  the  S  corner  from  that  in  either  the  NW  or  NE 
corner,  it  cannot  discriminate  between  the  NW  and  NE  regions  due  to  the  fact  they 
possess  identical  second-moment  properties.  As  a  partial  remedy  to  this  situ¬ 
ation  we  have  devised  a  discriminant  that  employs  both  correlation  and  edge 
density  information.  Since  this  discriminator  appears  sufficiently  interesting 
itself,  we  have  provided  details  on  its  implementation  in  the  Appendix.  Using 
this  correlation/edge  density  discriminator  some  degree  of  success  has  been 
achieved  in  discriminating  between  the  NW  and  NE  regions  as  illustrated  by 
the  results  in  Fig.  19d.  The  results  are,  however,  generally  inferior  to  the 
performance  of  the  log-likelihood  discriminator. 

In  order  to  assess  the  sensitivity  of  the  performance  of  the  log-likeli¬ 
hood  discriminator  to  our  modeling  assumptions,  we  have  applied  the  three  designs 

3h 


described  in  Table  1  to  the  original  image  in  Fig.  39c.  In  all  cases  we  have  taken 
d=3.  The  results  are  illustrated  in  Fig.  20  and  provide  some  indication 


of  how  performance  depends  upon  the  choice  of  filtering  function.  In  Fig.  20b 
we  illustrate  the  result  under  the  assumption  of  weak  interference  (Y^=6dB, 
p^-0. 0).  The  boundaries  between  the  three  different  texture  regions,  which 
are  actually  straight  lines,  are  irregular  and  small  patches  in  each  region 
have  been  misclassified.  In  Fig.  20c  we  illustrate  performance,  as  in  Fig.  19b, 
for  moderate  interference  (y^=0dB,  p^=0.0).  The  boundaries  between  regions 
are  smoother  and  the  misclassified  regions  are  smaller.  Finally,  in  Fig.  20d 
we  illustrate  performance  for  strong  interference  (y^=6dB,pi=0. 5).  The  bound¬ 
aries  are  now  much  straighter  and  the  misclassified  patches  have  disappeared. 
However,  the  point  of  intersection  of  the  three  boundaries  has  become  ill- 
defined  as  a  result  of  the  additional  smoothing  introduced  to  reduce  the 
interference.  These  results  are  useful  in  illustrating  the  tradeoffs  in 
choosing  the  strength  of  the  interference. 

VIII.  Summary  and  Conclusions: 

We  have  described  a  new  approach  to  texture  discrimination  which  appears 
to  offer  considerable  improvement  over  existing  approaches  under  specific, 
although  realistic,  stochastic  modeling  assumptions.  While  initial  results 
have  been  rather  encouraging,  much  more  work  remains  in  establishing  the 
efficacy  of  this  approach. 

For  example,  we  have  described  a  Wiener  filtering  approach  to  the  2-D 
digital  filter  implementation  of  the  class-conditional  generalized  log-like¬ 
lihood  functionals.  Clearly,  other  approaches  to  the  digital  filter  design 
problem  are  possible  within  the  general  structure  of  the  texture  discrim¬ 
inant  proposed  here.  Several  alternative  approaches  are  presently  under 
investigation. 
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A  more  fundamental  difficulty  with  this  approach,  however,  is  the  re¬ 
quirement  for  knowledge  of  the  model  ..rameters  for  each  of  the  texture 
classes  known  to  be  acting  in  any  discrimination  application.  This  know¬ 
ledge  is  required  in  constructing  the  inputs  g^( i , j ;d) ,k=0,l, . . . ,K-1,  to 
the  digital  filters  which  generate  the  class-conditional  log-likelihood  function¬ 
als.  We  assumed  this  knowledge,  for  example,  in  the  experiments  described  in 
the  preceding  section.  In  practice  it  may  be  possible  to  obtain  crude  esti¬ 
mates  of  these  parameters  either  on  the  basis  of  a  priori  knowledge  or  derived 
from  the  data  itself.  These  estimates  may  be  sufficiently  accurate  to  pro¬ 
vide  useful  discrimination  performance.  An  alternative  may  be  to  use  a 
range  of  "prototype"  stochastic  texture  models  which  span  the  range  of  textures 
of  interest.  Use  of  these  "prototype"  texture  models  in  constructing  corres¬ 
ponding  class -conditional  log-likelihood  functionals  may  result  in  useful  discrim¬ 
ination  ability.  Both  of  these  approaches  are  being  pursued. 

Finally,  the  proposed  texture  discriminant  is  intimately  related  to  our 
modeling  assumptions.  When  applied  to  realizations  of  stochastic  models  for 
which  it  was  developed,  the  performance  is  excellent.  Assessment  of  the  true 
value  of  this  approach,  however,  will  require  relative  performance  evaluation 
vis-a-vis  existing  approaches  on  real-world  data.  This  relative  evaluation 
should  be  respect  to  both  accuracy  and  computational  cost. 
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APPENDIX 


CORRELATION/EDGE  DENSITY  DISCRIMINATOR 

The  correlation/edge  density  discriminator  bases  its  decisions  on  maxi¬ 
mum  likelihood  estimates  p,A  of  the  correlation  of  points  separated  by  a  distance  d 
and  the  density  of  edges  respectively.  A  linear  discriminator  has  been  used  to  sepa¬ 
rate  the  p,X  plane  into  appropriate  regions.  This  discriminator  parallels 
the  log-likelihood  discriminator  to  the  extent  that  observations  of  pairs 
of  points  within  a  window  sure  assumed  independent  and  furthermore  the  window¬ 
ing  operation  is  replaced  by  a  weighted  window  implemented  as  a  2-D  recursive 
digital  filter.  Here,  however,  the  likelihood  functionals  are  maximized  over 
a  continuum  of  values  rather  than  a  finite  set  of  hypothesis. 

The  maximum  likelihood  estimate  of  the  correlation  will  be  formed  under 
the  assumption  that  observations  at  pairs  of  points  are  zero  mean  jointly 
Gaussian  random  variables.  The  p.d.f.  of  N  independent  observations  of  two 
jointly  Gaussian  random  variables  conditioned  upon  their  correlation  and 

4 

variance  is  given  by 
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Here  the  N  points  under  consideration  would  be  the  set  of  pairs  of  pixels 
within  a  window  W  and  separated  by  a  distance  d  while  the  random  variables 
x  and  y  correspond  to  a  pair  of  values  taken  on  by  each  pair  of  pixels. 

The  values  of  p  and  a  for  which  the  above  quantity  is  maximized  are  their 
maximum  likelihood  estimates.  It  will  prove  more  convenient  to  maximize  the 
logarithm  of  the  above  as  given  by 
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(A-2) 


l(p,o)®-N£n(2iro2)-  ^  in(l-p2) 


‘  “sVTi-S7) 


Upon  setting  the  partials  of  the  above  expression  with  respect  to  p  and  o  to 
zero  and  solving  for  p  and  a  we  obtain  the  estimates 


and 


•  N 

02=  l  (x2+y? )/2N  , 

i=l  1  1 

N 

P  =  l  x  y  /No2  . 


(A- 3a) 


(A- 3b) 


As  in  the  log-likelihood  discriminator  the  summations  above  over  a  window  were 
replaced  by  a  low  pass  filter. 

The  maximum  likelihood  estimate  of  the  edge  density  will  be  made  under  the 
assumption  that  the  occurrence  of  edges  along  an  arbitrarily  placed  line  can 
be  described  as  event  times  of  a  Poisson  process.  Then  the  probability  of 
observing  intervals  of  length  d  with  no  events  and  Ng  intervals  with  one 
or  iliore  events  under  the  assumption  that  the  edge  density  is  A  is  given 
by  [15] 

(N  +N  ) !  N 

P(Nq,nJa)  =  N°tNet  [exp(-Ad)]  0  [l-exp(-Ad)]  e  .  (A-U) 

0  e 

It  will  prove  convenient  to  work  with  the  logarithm  of  the  above  as 


(N  +N  )! 

MN0,NjA)  =  In  -N-U1-e1  ■  +  N  (-Ad)  +  Ne  £n[l-exp(-Ad) ]  (A-5) 

0  e 


Setting  the  derivative  of  the  above  with  respect  to  A  to  zero  and  solving  for 
A  yields 


X  =  d  £n  S0’ 


(A-6) 


where  8.  =  Nn/(N*+N  )  is  the  fraction  of  intervals  with  no  events.  In  fact, 
u  u  u  e 

the  quantity  used  in  the  discriminator  was  not  A  but  rather  an  estimate  of  0 


38 


formed  by  passing  the  output  of  a  simple  edge  detector  through  the  lorqmss 
filter  used  for  the  log -likelihood  discriminator.  The  edge  detector  simply 
applies  a  threshold  to  the  difference  between  the  gray  levels  at  pairs  of 
pixels  separated  by  distance  d.  Here  the  value  of  d  is  the  same  as  that  used 
by  the  log-likelihood  discriminator  (3  pixels)  while  the  lowpass  filter  is 
the  Viener  filter  designed  under  the  assumption  of  moderate  interference. 
Additional  details  can  be  found  in  [29]. 
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Figure  1 


Realizations  of  Two  Random  Fields  Possessing 
Identical  Power  Spectral  Densities; 

L,  X=0.08,  p=0.0;  R,  X=0.l6,  p=0.5 


Figure  2 

Selected  Realizations  of  Random  Fields 
Possessing  Identical  Edge  Density; 

W,  Rectangular  v=l,  X=0.08,  p=0.0, 

N,  Rectangular  v=°°,  X=0.08,  p=0.0, 
SE,  Polygonal  v=l,  X=0.32,  p=0.0. 


Rotation  of  Cartesian  Coordinate  Axes 


g.)  x=0.05,  p=-0.9  h.)  X=0.05,  P=0. 0  i. )  X=0.05, 


Figure  h 

Selected  Realizations  of  Random  Field 
Generated  by  Poisson  Partitions 


g. )  X»0.05,  p«-0.9  h.)  X«0.05,  p“0,0 


i.)  x=0.05,  p=0.5 


Figure  5 


Selected  Realizations  of  Random  Field 
Generated  by  Periodic  Partitions 


Figure  7 

of  Directed  Line  Segment 


H8 


Figure  8 

Selected  Realizations  of  Random  Field 
Generated  by  Polygonal  Partitions 


Figure  9 


Autocorrelation  Function  and  Power  Spectral  Density 
of  2-1)  Pander.  Checkerboard  Process  Generated  by 
Periodic  Partitions 


g. )  A  | | u | |=8.0,  o=-0. 9 


h. )  A  | |u| 1=8.0,  d=0.0 


i .  )  A 


1=8.0,  c=0.5 


Figure  10 

Selected  Joint  Probability  Density  Functions 
for  Rectangular  Partition  Process,  v=l 
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a-)  i Jul !=0.S.  p=-0.9 


u. )  i  iuj  |=l).D,  .=0.0 


|u|  '=0.5,  ,  =0.5 


d. )  \e  l  !u|  i  =  2.0.  .=-0.9  e.i  el'V,  •  =o.O  f.)  »e!|u||=2.0,  e=0.5 


Figure  11 

Intensity  Plots  of  Logarithm  of  Selected 
Joint  Portability  Density  Functions 
for  Rectangular  Partition  Process,  v=l 
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»•)  Ael lul 1-0.5.  p*-0. 9 


M  1*0-5.  P'0.0 


X 


d.)  *_l |u| 1=2.0.  p=-0.9 


c.)  »e||u| 1-0.5.  P*0.5 


0 


f.)  X.llul 1-2.0,  p-o.5 


Figure  13 

Intensity  Plots  of  Logarithm  of  Selected  Joint 
Probability  Density  Functions  for  Polygonal 
Partition  Process 
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ypical  2-D  Random  Field  Models  of  Texture 
Rectangular  Partition  Process  with  A=0.l6 
and  Various  Values  of  p  and  v 
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Figure  16 

Point  Spread  Function  of  Composite  Filter 


NORMALIZED  RESPONSE,  IN  dB 


RADIAL  FREQUENCY  0/2t  RADIAL  FREOUENCY  Ii/2* 

a.)  Weak  Interference  Model  b.)  Moderate  Interference  Model 


RADIAL  FREQUENCY  Q/Zw 


c.)  Strong  Interference  Model 


Figure  17 

Illustration  of  'typical  Pover  Spectral  Densities 
And  Associated  Wiener  Filter  Response 
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)  Dp oired 
Strong 


Response  for  f.  )  Actual  Response  for 

Interference  Strong  Interference 


Figure  l8 


Frequency  Responses  of  Desired  and  Actual  Filters 


a. )  Original 


b.)  Weak  Interference 


)  Moderate  Interference  d. )  Strong  Interference 


Figure  20 

Performance  of  Log-Likelihood  Texture  Discriminator 
for  Various  Assumed  Interference  Levels 


